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We have used our new technique for fully numerical evolutions of orbiting black-hole binaries 
without excision to model the last orbit and merger of an equal-mass black-hole system. We track the 
trajectories of the individual apparent horizons and find that the binary completed approximately 
one and a third orbits before forming a common horizon. Upon calculating the complete gravitational 
radiation waveform, horizon mass, and spin, we find that the binary radiated 3.2% of its mass and 
24% of its angular momentum. The early part of the waveform, after a relatively short initial 
burst of spurious radiation, is oscillatory with increasing amplitude and frequency, as expected 
from orbital motion. The waveform then transitions to a typical 'plunge' waveform; i.e. a rapid 
rise in amplitude followed by quasinormal ringing. The plunge part of the waveform is remarkably 
similar to the waveform from the previously studied 'ISCO' configuration. We anticipate that the 
plunge waveform, when starting from quasicircular orbits, has a generic shape that is essentially 
independent of the initial separation of the binary. 

PACS numbers: 04.25.Dm, 04.25. Nx, 04.30.Db, 04.70.Bw 



I. INTRODUCTION 

The study of the late orbital stage of black-hole bi- 
naries is of particular interest because they are thought 
to be the most likely sources detected by gravitational 
wave observatories such as LISA |l( and LIGO 2] (which 
is now reaching its design sensitivity). Besides, the two 
body problem in General Relativity is, in itself, a gen- 
uinely interesting theoretical problem. Studies of fami- 
lies of binary-black-hole initial data in quasicircular or- 
bits set the periods of the innermost stable circular or- 
bit (ISCO) to 37M for the Bowen-York family of initial 
data 0, 3, 57 M for the Thin-Sandwich family H @, 
and 49-M from third-order post-Newtonian approxima- 
tion (3PN) 0. Evolution of binary black holes from these 
locations lead to plunge motion 0, 0, ^j; performing 
a fraction of an orbit before a common horizon encom- 
passes the two black holes (i.e. forming a single, distorted 
black holes). 

There have been several remarkable advancements in 
Numerical Relativity in the past few years, and, particu- 
larly since the work of 0, Il2l ^| , it now seems possible 
to evolve orbiting black-hole binaries out from arbitrary 
distances to the merger and ringdown. The first major 
breakthrough in numerical evolutions of these systems 
was reported by Briigmann et. al. |14|. Using a careful 
choice of corotating shift and singularity excision, they 
were able to evolve a black-hole binary, starting from 
initial data for a quasicircular binary, for more than a 
complete orbit. Their work was recently verified by Di- 
ener et. al. ^{|. However, in both cases it was not pos- 
sible to extract the waveform from the merger. The first 
fully numerical evaluation of the waveform from an or- 
biting black-hole binary was reported by Pretorius [l3j . 
Pretorius evolved a system in which two scalar fields col- 
lapsed to form individual black holes, which then formed 



a merging elliptical binary. The evolution used a direct 
discretization of a second order in time system with sin- 
gularity excision. Recently, a new technique 0, 0] has 
been developed for evolving black-hole binaries that uses 
the more conventional BSSN [lti, HtI Il8l| system of equa- 
tions (which are first-order in time) . This new technique 
is based on the puncture approach, but allows the sin- 
gular punctures to move across the grid. Singularity ex- 
cision is not required and the new system has the ad- 
vantages that it is very easy to implement and appears 
to be very accurate. This is explicitly demonstrated in 
Ref. where this technique was applied to study the 
merger kicks of unequal-mass black hole binaries verify- 
ing those previously computed in Ref. [2fj| . 

We use this new 'moving puncture' approach to evolve 
the same initial configuration as in |14| and confirm that 
the system does indeed undergo more than a full orbit 
before a common horizon forms. We were also able to 
accurately extract the waveform and final horizon pa- 
rameters. We find very good agreement in the radiated 
energy, angular momentum, and merger time between 
those calculated from the horizon properties and from 
the waveform. 



II. FORMULATION 

Our simulations of orbiting black-hole binaries are 
based on a modification to the standard puncture ap- 
proach. In the puncture approach |2l| the metric on the 
initial slice is given by 22] j a b — {ipBL + u) 4 5 a b, where 
ipBL = 1 + J2i=i m i/{^ r i) is the Brill-Lindquist confor- 
mal factor, is the mass parameter of puncture i, ri is 
the coordinate distance to puncture i, and u is finite on 
the punctures. 

In the standard puncture approach the locations of the 
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punctures are fixed (one imposes that the shift vanishes 
at the puncture location) , and the singular behavior (i.e. 
i/jbl) in the metric is handled analytically. A conse- 
quence of fixing the punctures is that the coordinates 
become highly distorted, and this, in turn, causes nu- 
merical runs of orbiting black holes (without excision or 
corotation) to crash relatively early. In a recent paper we 
introduced a new technique for evolutions with moving 
punctures 

[12 (see [12 for an alternative implementa- 
tion). This new technique, which is based on the BSSN 
formulation of General Relativity [T^ . ll7llT8l |. does not re- 
quire either excision or a corotating shift. Our technique 
replaces the BSSN conformal exponent </>, which is infi- 
nite on the punctures, with the C 4 field \ — exp(— 40). 
This new variable, along with the other BSSN variables, 
will remain finite provided that one uses a suitable choice 
for the gauge. 

We obtained accurate, convergent waveforms by evolv- 
ing this system in conjunction with a modified 1 +log 
lapse, a modified Gamma-driver shift condition 0, |23| . 
and an initial lapse set to a — V'bl- The lapse and shift 
are evolved with 

d a 
d t p a 



and the 



-2aK, (1) 
B a , d t B a = 3/4<9 t f a - r)B a . (2) 



We use the LazEv framework [24( to numerically evolve 
this new system. Unlike in the fixed puncture approach, 
we do not reduce the order of finite differencing near the 
punctures. We use the standard centered fourth-order 
stencils for all derivatives except for the advection terms 
(i.e. terms of the form f3 l di) where we use upwinded 
fourth-order stencils. These stencils were modified near 
the boundary. At the second point from the boundary 
we use fourth-order centered stencils for all derivatives, 
and at the first point from the boundary, we use second- 
order centered stencils. We use the standard fourth-order 
Runge-Kutta algorithm for the time evolution and radia- 
tive boundary conditions for all evolved variables. 



III. INITIAL CONFIGURATION 

Following [12 we choose black hole initial data from a 
quasicircular sequence |25| with parameters 

m/M = 0.47656, P/M = 0.13808, Y/M = ±3.0, 
L/M = 9.32, J/M 2 = 0.82843, Mfl = 0.054988, 

where m is the mass of each single black hole, M = 1 
provides the scale, P is the magnitude of the linear mo- 
menta (equal but opposite and perpendicular to the line 
connecting the holes), (0, Y, 0) is the coordinate location 
of the punctures, L is the proper distance between the ap- 
parent horizons, J is the total angular momentum, and f2 
is the orbital frequency. We use the Brandt-Brugmann 
approach along with the BAM_Elliptic [22|, |2g Cactus 
thorn to solve for these initial data. The total ADM 
mass for this configuration is Madm — 0.98461Af. 



The initial choice for the lapse is a - 
initial choice for the shift is (5 % = B l = 0. 

We evolved these data with grid resolutions of M/21, 
M /24, and M /27; and gridsizes of 224 2 x 448, 256 2 x 512, 
and 288 2 x 576 respectively (we exploited the 7r-rotation 
symmetry about the z axis and reflection symmetry 
about the equatorial plane to halve the number of grid- 
points in the x and z directions). We used a 'multiple 
transition' Fisheye transformation, which is an exten- 
sion of the 'transition' Fisheye transformation [23l |27| , 
to place the boundaries at 114M. The 'multiple transi- 
tion' transformation has the form R = C r, where R is 
the physical radius corresponding to the coordinate ra- 
dius r, 



C 



E, . i cosh((r + r0i)/si 
K t /r\og 

i=\ 



cosh((r - r0i)/si 



(oj_i - a l )s i 
2 tanh(r0i/sj 



(3) 
(4) 



n is the number of transitions, a; is the deresolution pa- 
rameter in region i, ao is the central resolution, rOi is the 
center of the ith transition, and Si is the width of the 
ith transition. For these runs we used the parameters, 
n = 2, a = 1, cii = 5, a 2 = 30, rOi = 5, r0 2 = 7.5, 
si = S2 = 0.75. We also evolved these data with a reso- 
lution of M/21 and a gridsize of 288 2 x 576, which placed 
the boundary at 176. 6M, to quantify the dependence of 
the results on the location of the boundary. For this run 
we used the Fisheye parameters, n = 1, ao = 1, a± = 25, 
rOi = 7.0, si = 25. 



IV. RESULTS 

We used Jonathan Thornburg's AHFinderDirect 
thorn 28] to find apparent horizons. The (coordinate) 
time of first appearance of the common horizon depends 
on resolution, with T ca h = 111M for the M/21 run, 



113M (estimated) for M/24, and T a 



114.3 



(estimated) for the M/27 run. An extrapolation of 
these data to infinite resolution puts the appearance 
of the first common horizon at T ca h — 125AZ. The 
common horizon has an irreducible mass of Mi rr — 
(0.8848 ±0.0002)M and specific spin (measured from the 
distortions of the horizon) of a — J-^/M^ — 0.688 ± 
0.001. These parameters correspond to a horizon mass 
of M H = (0.952 ± 0.002)M and angular momentum 
J H = (0.6232 ± 0.003)M 2 . Hence (3.3 ± 0.2)% of the 
mass and (24.7 ± 0.4)% of the angular momentum are 
radiated away. 

Figure^shows the tracks of the punctures, the individ- 
ual horizons every 10M of evolution, and the first com- 
mon horizon. The plot was generated using the M/21 
resolution run with the boundaries at 176. 6M. In ad- 
dition the plot shows the puncture trajectory for the 
M/27 resolution run. Note that the binary completes 
one and a third orbits before the common horizon forms 
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FIG. 1: The puncture trajectories and apparent horizon pro- 
files on the xy plane for the M/21 run. The solid and dotted 
spirals are the puncture trajectories, the solid and dotted el- 
lipsoids are the individual apparent horizons (at every 10M 
of evolution), the dot-dash 'peanut shaped' figure is the first 
detected common horizon, and the dashed spiral is the punc- 
ture trajectory for the M/27 run. The initial growth of the 
individual apparent horizon is due to the non-ideal (vanish- 
ing) initial data for the shift. Note that we track the puncture 
positions throughout the evolution. The period of the last or- 
bit is around 62M. The last orbit begins when the punctures 
are located at 2.6M from the origin. 



(although we caution the reader that tracks are gauge de- 
pendent). The period of the last orbit is approximately 
62M. The puncture trajectories were calculated by in- 
tegrating d t x punct = -^ punct , where f} punct is the inter- 
polated value of the shift on the puncture (the puncture 
never lies on a gridpoint). The last orbit begins when 
the punctures are located at 2.6M from the origin in the 
coordinate conformal space. 

We use the Zorro thorn 0> H3 to calculate tpi an d 
decompose it into spherical harmonics of spin weight —2. 
The two dominant modes are (£ — 2,m — +2) and (£ = 
2,m = —2), where the coefficient of the two modes are 
complex conjugates. Figure [3 shows the (£ = 2, m = 2) 
mode of ip4 at r = 20M for the three resolutions and 
a convergence plot of these data. The waveforms con- 
verge to fourth-order up to t — 141M (the convergence 
rate past t = 141M is obscured by the large phase error, 
but is better than second-order). The oscillation in -04 
around t = 26M is due to spurious radiation in the ini- 
tial data. This spurious radiation quickly leaves the sys- 
tem and becomes smaller than the orbital waveform at 
t ~ 50M . Hence the radiation from the last orbit, which 
begins at t ~ 50M as seen by our observer located at the 
radial coordinate r = 20M (see Figure^, is essentially 
uncontaminated. We find that the (£ = 2, m = 2) quasi- 
normal mode, for the M/21 resolution run with distant 



FIG. 2: The (£ = 2, m = 2) mode of ^4 at r = 20M. The 
top plot shows the waveforms for central resolutions of M/21, 
M/24, and M/27. The bottom plot shows the differences 
4> 4 (M/2l) - V4(M/24) and 4> 4 {M/24) - 4> 4 (M/27), with the 
latter rescaled by 1.879 to demonstrate fourth-order conver- 
gence. Note the spurious radiation around t = 26M. 



boundaries, has a frequency of uj/a.20 = 0.549 ±0.001 
("20 = 0.954 is the average value of the lapse at r = 20M 
at late times). The reported error is from the fit to a 
damped sinusoidal function and does not include finite 
difference errors. This frequency corresponds |29j to a 
specific spin of a = Jn/M^ = 0.673 ±0.002. We can also 
estimate when the first common horizon forms by calcu- 
lating the offset between the tallest peak in the plunge 
waveform for these data and that of the 'ISCO' (as deter- 
mined by the effective potential method) initial data. In 
a previous paper ^l| we found that a common horizon 
forms in the 'ISCO' case at t = 19.333M. From these 
offsets we estimate that the first common horizon should 
form at 111.5M for the M/21 run, 112. 9M for the M/24 
run, and 113.7M for the M/27 run. These numbers are 
within 0.6M agreement with those determined directly 
from the apparent horizon finder. An extrapolation of 
these estimates for the formation of the common horizon 
yields T cah = 121.5M . 

We calculated the radiated energy and angular mo- 
mentum from ipi at r = 15M, r — 20M, r = 25M, and 
r = 30M. We then extrapolated these data (based on 
a least squares fit versus p = 1/r) and found that ex- 
trapolated radiated energy and angular momentum were 
(3.18±0.2)% and (24.3±2)% respectively. Notably, these 
results are in excellent agreement with those calculated 
from the horizon mass and angular momentum. Tabled 
summarizes the physical parameters extracted indepen- 
dently from the horizon and waveform. 

The plunge part of the waveform shows remarkable 
similarities with the waveform from the ISCO. Figure |3| 
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FIG. 3: The real part of the [i = 2, m = 2) mode of the tpi 
at r = 20M from this 'last orbit' configuration and from the 
ISCO configuration. Note the near perfect overlap once the 
ISCO waveform has been translated by At/Af = 95 
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FIG. 4: The (£ = 2,m = 2) component of the strain. Both the 
+ and x mode are shown. The early time strain is dominated 
by spurious radiation (from the initial data) up to t = 55M. 
Afterwords, the strain shows a gradual transition from orbital 
motion to a plunge waveform. This transition is less distinct 
that that in ip4. 



shows the real part of the (£ = 2, m = 2) mode for these 
two configurations. Note that, after a translation, there 
is near perfect overlap of the late-time waveform. 

The gravitational strain h is related to "04 by ip^ = 
-l/2d?h(t). In Fig.Hwe show the {I = 2, m = 2) com- 
ponent of both polarizations of the strain. The early 
part of the strain (t < 55M) is dominated by the spu- 
rious radiation of the initial data. Note that the strain 
amplitude and frequency shows a far more gradual tran- 
sition from an orbital inspiral type waveform to a plunge 
type waveform than ip4 and seems better suited to match 
to post-Newtonian waveforms. 
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FIG. 5: The Hamiltonian constraint violation at t = 70 M 
along the y-axis for the M/24 and M/27 runs (the latter 
rescaled by (27/24) 4 . The punctures crossed the y-axis for 
the second time at t = 64M. Note the near perfect fourth- 
order convergence. Points contaminated by boundary errors 
have been excluded from the plot. The high frequency vio- 
lations near the numerical coordinate y = ±8 are due to the 
extreme fisheye deresolution near the boundary, and converge 
with resolution. 



TABLE I: Results of the evolution as determined from the 
waveform and the remnant horizon. 

Method Erad/MADM Jrad/JADM T cah / M 0, / ' M H 

Radiation (3.18 ±0.2)% (24.3 ±2)% - 121 0.673 ± 0.002 
Horizon (3.3 ± 0.2)% (24.7 ± 0.4)% w 125 0.688 ± 0.001 



To demonstrate consistency with the General Relativ- 
ity field equations, we calculated the Hamiltonian con- 
strain violation. The constraint converges to fourth-order 
outside of a small region surrounding the puncture (the 
Hamiltonian constraint violation on the nearest neigh- 
boring points to the puncture are roughly independent of 
resolution, but this non-converging error does not prop- 
agate outside the individual horizons). Figure [S] shows 
the Hamiltonian constraint violation along the y-axis at 
t = 70M (about AM after the punctures cross the y 
axis for the second time). The constraint is fourth or- 
der convergent everywhere but at points contaminated 
by boundary errors (these points have been removed from 
the plot). 



V. DISCUSSION 

Using our new technique that allows punctures to move 
in the numerical grid, we have succeeded in accurately 
computing the last orbit of a black-hole binary, obtained 
waveforms and extracted relevant physical information 
such as energy and angular momentum radiated, appar- 
ent horizon geometry, and orbital parameters. These re- 
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suits are consistent with each other, as summarized in 
Tabled We also note the interesting fact that the plunge 
part of the waveforms (corresponding to the highest am- 
plitude region in ip^) is roughly insensitive to the initial 
separation of the holes when starting from a quasicircular 
orbit (see Fig. Inhere and Fig. 30 in Ref. Q) 

When evolving the same initial configuration, 
Briigmann et. al. did not find a common horizon 
(they evolved to 185Af), and concluded that the binary 
must have undergone at least one orbit since the orbital 
period of the initial configuration is around 120M. 
Recently, Diener et. al. [15| evolved this system with 
various choices for the gauge (all containing a corotating 
shift) and concluded that a common horizon forms at 
about 120M - 125M (after the binary completes a full 
orbit). They get these estimates by evolving with the 
extremely high resolutions of M/66 and M/80, and 
extrapolating to the continuum limit. They compute 
the trajectory of the apparent horizons and find that 
the orbital period of the last full orbit is 59 M. In 
both 01 an d waveform extraction was not possible. 
In our work we find, like Diener et. al., that the common 
horizon forms between 120M and 125M, and that the 
orbital period of the last orbit is approximately 62M. 
However, we only required a resolution as high as M/27, 
and could calculate accurate waveforms. 



Aside from the duration of the last orbit, we estimate 
that the initial separation of the black holes in the final 
orbit is 5.2M in conformal coordinates. It is interesting 
to compare this separation with those of several ISCO de- 
terminations. For Bowen-York initial data (as was used 
in this paper) the ISCO separation 0, 0| is 2.3M, for 
Thin-Sandwich data 0, E| it is 3.25M, and an estimate 
of the 3PN @ ISCO puts the separation at 4.24M . Ob- 
viously, radiation reaction, which was not taken into ac- 
count in those computations, leads to radial motion that 
plays an important role in the dynamics of the last orbit 
of black-hole binaries. 
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